Gene Expression Analysis

Fetal colon cell line FHC | Necroptosis, Ferroptosis & Pyroptosis

I. Preliminaries

Loading libraries

library("tidyverse")
library("tibble")
library("msigdbr")
library("ggplot2")
library("ensembldb")
library("purrr")
library("magrittr")
library("matrixStats")
library("dplyr")
library("grex")
library("gplots")
library("RColorBrewer")

Constants

DATA_DIR <- "data/FHC/"

Loading the Expression Data

The expression data are taken from this study: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE232211

Download (1) the RNA-seq normalized counts matrix and (2) the human gene annotation table by running:

wget https://www.ncbi.nlm.nih.gov/geo/download/?type=rnaseq_counts&acc=GSE232211&format=file&file=GSE232211_norm_counts_TPM_GRCh38.p13_NCBI.tsv.gz -P data/FHC/
wget https://www.ncbi.nlm.nih.gov/geo/download/?format=file&type=rnaseq_counts&file=Human.GRCh38.p13.annot.tsv.gz -P data/FHC/

fhc.expression <- read.delim(paste0(DATA_DIR, "GSE232211_norm_counts_TPM_GRCh38.p13_NCBI.tsv"), as.is = TRUE, header = TRUE, row.names = 1)
fhc.expression <- rownames_to_column(fhc.expression, "gene_id")
fhc.expression

mapping_table <- read.delim(paste0(DATA_DIR, "Human.GRCh38.p13.annot.tsv"), as.is = TRUE, header = TRUE, row.names = 1)
mapping_table <- rownames_to_column(mapping_table, "gene_id")
mapping_table <- subset(mapping_table, EnsemblGeneID != "")
mapping_table

Since the gene set for the RCD-related regulators (taken from RCDdb) uses Ensembl Gene IDs, we need to map the accessions.

fhc.expression.mapped <-  right_join(fhc.expression, mapping_table, by = join_by(gene_id == gene_id)) %>% distinct(EnsemblGeneID, .keep_all = TRUE)
rownames(fhc.expression.mapped) <- fhc.expression.mapped$EnsemblGeneID
fhc.expression.mapped = fhc.expression.mapped[,!(names(fhc.expression.mapped) %in% colnames(mapping_table))]
fhc.expression.mapped <- rownames_to_column(fhc.expression.mapped, "gene_id")
fhc.expression.mapped

Exploratory Data Analysis

We load the gene sets from RCDdb: https://pubmed.ncbi.nlm.nih.gov/39257527/

RCDdb <- "data/RCDdb/"

Necroptosis

Load the gene set.

genes <- read.csv(paste0(RCDdb, "Necroptosis.csv"))
genes$gene_id <- cleanid(genes$gene_id)
genes <- distinct(genes, gene_id, .keep_all = TRUE)
genes <- subset(genes, gene_id != "")
genes

Get the TPM for the genes in the gene set.

tpm.df <- fhc.expression.mapped %>% dplyr::filter(gene_id %in% genes$gene_id)
tpm.df <- left_join(tpm.df, genes %>% dplyr::select(gene_id, gene), by = c("gene_id" = "gene_id"))
rownames(tpm.df) <- tpm.df$gene
tpm.df <- tpm.df %>% dplyr::select(c(1:4))
tpm.df <- subset(tpm.df, select = -c(gene_id) )
tpm.df <- tpm.df[ order(row.names(tpm.df)), ]
tpm.df

Plot the results.

tpm.matrix <- as.matrix(tpm.df)
heatmap.2(tpm.matrix, srtCol=360, cellnote = tpm.matrix, dendrogram="none", Colv=FALSE, Rowv=FALSE,
          col=brewer.pal(n = 9, name = "BuPu")[5:9], trace="none", key = FALSE, lwid=c(0.1,4), lhei=c(0.1,4),
          cexCol=1, cexRow=0.75, symm = TRUE)

Ferroptosis

Load the gene set.

genes <- read.csv(paste0(RCDdb, "Ferroptosis.csv"))
genes$gene_id <- cleanid(genes$gene_id)
genes <- distinct(genes, gene_id, .keep_all = TRUE)
genes <- subset(genes, gene_id != "")
genes

Get the TPM for the genes in the gene set.

tpm.df <- fhc.expression.mapped %>% dplyr::filter(gene_id %in% genes$gene_id)
tpm.df <- left_join(tpm.df, genes %>% dplyr::select(gene_id, gene), by = c("gene_id" = "gene_id"))
rownames(tpm.df) <- tpm.df$gene
tpm.df <- tpm.df %>% dplyr::select(c(1:4))
tpm.df <- subset(tpm.df, select = -c(gene_id) )
tpm.df <- tpm.df[ order(row.names(tpm.df)), ]
tpm.df

Plot the results.

tpm.matrix <- as.matrix(tpm.df)
heatmap.2(tpm.matrix, srtCol=360, cellnote = tpm.matrix, dendrogram="none", Colv=FALSE, Rowv=FALSE,
          col=brewer.pal(n = 9, name = "BuPu")[5:9], trace="none", key = FALSE, lwid=c(0.1,4), lhei=c(0.1,4),
          cexCol=1, cexRow=0.75, symm = TRUE)

Pyroptosis

Load the gene set.

genes <- read.csv(paste0(RCDdb, "Pyroptosis.csv"))
genes$gene_id <- cleanid(genes$gene_id)
genes <- distinct(genes, gene_id, .keep_all = TRUE)
genes <- subset(genes, gene_id != "")
genes

Get the TPM for the genes in the gene set.

tpm.df <- fhc.expression.mapped %>% dplyr::filter(gene_id %in% genes$gene_id)
tpm.df <- left_join(tpm.df, genes %>% dplyr::select(gene_id, gene), by = c("gene_id" = "gene_id"))
rownames(tpm.df) <- tpm.df$gene
tpm.df <- tpm.df %>% dplyr::select(c(1:4))
tpm.df <- subset(tpm.df, select = -c(gene_id) )
tpm.df <- tpm.df[ order(row.names(tpm.df)), ]
tpm.df

Plot the results.

tpm.matrix <- as.matrix(tpm.df)
heatmap.2(tpm.matrix, srtCol=360, cellnote = tpm.matrix, dendrogram="none", Colv=FALSE, Rowv=FALSE,
          col=brewer.pal(n = 9, name = "BuPu")[5:9], trace="none", key = FALSE, lwid=c(0.1,4), lhei=c(0.1,4),
          cexCol=1, cexRow=0.75, symm = TRUE)

---

1. De La Salle University, Manila, Philippines, gonzales.markedward@gmail.com↩︎

2. De La Salle University, Manila, Philippines, anish.shrestha@dlsu.edu.ph↩︎